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Abstract 

Let AF be the free energy difference between two equilibrium states of a sys- 
tem. An established method of numerically computing AF involves a single, 
long "switching simulation" , during which the system is driven reversibly from 
one state to the other (slow growth, or adiabatic switching). Here we study a 
method of obtaining the same result from numerous independent, irreversible 
simulations of much shorter duration (fast growth). We illustrate the fast 
growth method, computing the excess chemical potential of a Lennard-Jones 
fluid as a test case, and we examine the performance of fast growth as a 
practical computational tool. 
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INTRODUCTION 



The development of numerical methods of computing the free energy difference AF = 
F B — F A between two equilibrium states of a system, A and B, is a field of active research0@, 
of significant practical importance in chemistry, physics, and biology. One widely applicable 
method - often called slow growth or adiabatic switching - utilizes a dynamical simulation 
during which the system evolves, maintained at a constant temperature, as an external pa- 
rameter A is slowly changed so as to drive the system from the initial to the final equilibrium 
state, A — > B. The evolution can occur as a sequence of discrete Monte Carlo steps, or as a 
time-continuous molecular dynamics trajectory. The slow growth method is often presented 
as a variant of thermodynamic integration^, with AF expressed as an integral over a contin- 
uous sequence of equilibrium states. However, slow growth can equally well be understood 
as the in silico exploitation of a basic result of statistical mechanics^: 

AF = W 001 (1) 

which states that the external work performed on a system during a reversible, isothermal 
process connecting two equilibrium states, is equal to the free energy difference between 
them. The slow growth method amounts to simulating such a "switching process" numeri- 
cally, then equating the work performed with the desired free energy difference. 

The subscript on emphasizes that the process carrying the system from A to B 
is reversible, hence (in principle) of infinite duration. During such a process, the system 
evolves through a continuum of equilibrium states connecting A to B. In practice, exact 
reversibility is never achieved, and one uses a finite-time simulation to estimate the free 
energy difference: 

AF « W T , (2) 

where W T denotes the work performed during a simulation of duration r. The accuracy of 
a slow growth computation depends on how close one is to the reversible limit; for short 
switching times r, the system can be driven far from equilibrium, resulting in a poor estimate 
of AF. 

In recent years, the following non- equilibrium work relation has been derived!'!, suggest- 
ing an alternative prescription for computing AF: 



exp(-(3AF) = exp(-/W r ). (3) 

This result pertains to an ensemble of simulations of a switching process of arbitrary duration 
r. The simulations are carried out independently, with initial conditions for the system 
sampled randomly from the canonical ensemble corresponding to the equilibrium state A. 
The work W T is computed for each simulation, and the overbar on the right denotes an 
average over the ensemble. The factor j3 is the inverse temperature of the initial equilibrium 
state, and also of the thermostat (see Section | below). Eq.|] can be rewritten as: 



AF= lim W*' , (4) 

N— >oo 



where 
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1 N 

-£exp(-/W r , n ) 

n=l 



(5) 



Here, W r n is the work performed during the n'th of iV independent simulations. We will 
refer to W*' N as the exponential average of W T over the iV simulations (in contrast with the 
ordinary average, which we will later introduce as W*' N ). We stress that the non-equilibrium 
work relation is valid for arbitrary r. This is perhaps surprising: even for short switching 
times - for which the system is driven far from equilibrium, and the work performed during a 
typical simulation gives a sorry estimate of the free energy difference - we can still compute 
AF from the set of work values {W T) i, W T> 2, ■ ■ • , W T>n , ■ ■ •} obtained from an ensemble of 
simulations. 

Since W* ,N —> AF as N — > oo (with r fixed), we are led to the following large- N 
approximation: 

AF « W? N . (6) 

Thus, after running many simulations of a given switching process, we can estimate AF 
as the exponential average of the values of work performed. We will refer to this proposed 
method of computing free energy differences as the fast growth method, to emphasize that 
it remains valid even for short switching times r. 

The aim of this paper is a study of the fast growth method. We have performed numerous 
switching simulations of the insertion of a single particle into a Lennard- Jones fluid, at 
various switching times r. We will use the results of these simulations both to present an 
illustration of principle, and to examine the utility of fast growth as a practical computational 
method. 

This paper is organized as follows. In Section | we briefly review the theoretical founda- 
tion of the fast growth method. In Section |TJ we introduce our test case - the evaluation of 
the excess chemical potential of a Lennard- Jones fluid - and we present details of the com- 
putation. In Section |TTT| we use our numerical results to illustrate the validity of fast growth. 
In Section |V| we use the same results to compare the slow and fast growth methods, given 
a fixed amount of computer time. Finally, in Section |V] we comment on the relationship 
between slow growth, fast growth, and a free energy formula based on the linear response 
approximation. 



I. THEORETICAL BASIS OF FAST GROWTH 

Here we introduce notation and briefly review the theoretical basis of the fast growth 
method. The emphasis will be on a presentation of the non-equilibrium work relation (Eq.^J) 
oriented specifically toward the numerical problem of computing free energy differences. For 
a more physical presentation, see Reffi 

Consider a system with a finite number of degrees of freedom, D. Let z denote a point 
in the 2.D-dimensional phase space of the system. Next, let A denote an externally con- 
trolled parameter (for instance, an applied field), and let H\(z) be the parameter-dependent 
Hamiltonian which gives the internal energy of the system, in terms of its microstate z. 
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Now consider two equilibrium states of the system, A and B, denned at the same tem- 
perature, T, but corresponding to different values of the external parameter, X A and X B . 
The partition function associated with state A is given by the integral (assumed finite) 

Z A = J dzexp[-(3H A {z)], (7) 

where (3 = (/c^T) -1 , and H A = H x a; and similarly for state B. The free energy difference 
between the two states is then 

AF = F B — F A = — In (8) 

It is this quantity that we want to compute. 

Next, suppose we choose a thermo stating scheme, a numerical prescription for generat- 
ing evolution under which our system samples the canonical ensemble associated with any 
equilibrium state (A, T). Intuitively, we can view the thermostating scheme as a way of 
"mocking up" the evolution of a system in contact with a heat reservoir. This evolution 
might take the form of a continuous-time trajectory, or a discrete sequence of Monte Carlo 
steps. Examples of such schemes include the Metropolis Monte Carlo algorithmic, Langevin 
dynamics, the Andersen thermostat!, and the Nose-Hoover thermostatic. Note that the 
thermostating scheme is itself parametrized by (A, T) (e.g. the Metropolis acceptance prob- 
ability takes the form e~ f3Hx ), and might be explicitly stochastic, relying on the generation 
of random numbers to produce the evolution. 

Imagine that we now run a simulation, first sampling the initial microstate of our system 
from the canonical ensemble corresponding to state A, then letting the system evolve in time 
under the thermostating scheme, as the value of A is externally changed - or "switched" - 
from \ A to X B , and the temperature of the thermostat is held fixed at T. For continuous- 
time evolution, let r denote the duration of the simulation; let X(t) specify the (externally 
imposed) parametric time-dependence, from A(0) = \ A to A(r) = \ B ; and consider the 
integral 

W T = I dtX(t)-^(z(t)), (9) 

where A = dX/dt, and z(t) represents the evolution of the system. Alternatively, if the 
evolution is implemented with an explicitly discrete algorithm (e.g. Metropolis Monte Carlo), 
then the corresponding quantity is 

W T = J2[H Xt+1 (z t )-H Xt (z t )}, (10) 
t=o 

where z t is the microstate of the system and X t the value of the parameter, after t steps; and 
t is the total number of steps in the simulation. In what follows, we will use the notation of 
continuous-time evolution, but the statements apply to explicitly discrete evolution as well. 

The quantity W T can be interpreted as the external work performed on the system over 
the course of the simulation.BMJrE-3 We will use this interpretation, though we warn the 
reader that Eq|] (or Eq.|Tl]) is by no means an accepted, textbook definition of work, and 
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the interpretation of W T as work may be subject to dispute. The results presented, however, 
do not depend on the validity of this interpretation. 

The trajectory z(t) which emerges from our simulation (and which in turn determines the 
value of W T ) depends on a string of random numbers: those involved in sampling the initial 
conditions, as well as those, if any, needed to generate the evolution itself. By repeating 
the simulation with different sets of random numbers, we obtain different realizations of 
the process we are simulating. We view these as representing independent samples from 
an infinite ensemble of all possible realizations.^ Eq.|3] then states that the average value 
of e~@ WT over the ensemble of possible realizations is equal to e~ l3AF , with AF as given by 
Eq||. In other words, e~^ Wr is an unbiased estimator of e~ l3AF . By collecting sufficiently 
many samples from the ensemble of realizations (i.e. by repeatedly running simulations) we 
can in principle compute the average of e~^ Wr to the desired degree of accuracy. This is the 
essence of the fast growth method. 

For a variety of thermostating schemes, Eq.|3| has been established as an exact result, 
valid for any finite switching time r. Precise conditions under which this result is rigorously 
true are presented in RefsHS. In particular, all four thermostating schemes mentioned 
earlier - Metropolis, Langevin, Andersen, Nose-Hoover - satisfy these conditions. 

Note that we do not automatically gain something for nothing with the fast growth 
method: while it allows us to compute AF without the bother of running a very long 
simulation, the price paid is that many short simulations might be needed. 

Finally, we mention that, in the limit r — > 0, Eq.|3] reduces to the well-known free energy 
perturbation identity^: 

exp(-/3AF) = (exp(-/3AH)) A , (11) 

where AH = H B — H A is the difference between the two Hamiltonians, evaluated at a given 
microstate; and (• • -)a denotes an average over microstates sampled from the canonical 
ensemble A. 



II. TEST CASE: EXCESS CHEMICAL POTENTIAL OF A LENNARD-JONES 

FLUID 

To illustrate the fast growth method, we chose to compute the excess chemical potential 
of a three-dimensional Lennard- Jones fluid. This is equal to the free energy change associated 
with "turning on" the interactions between a single particle and the remainder of the fluid, 
at constant volume and temperatureB. This process can be viewed as the transfer of one 
particle from an ideal gas to the Lennard- Jones fluid; the excess chemical potential is defined 
as the chemical potential of the fluid relative to that of the ideal gas. Following the notation 
of Mon and Griffiths!! let Um denote the potential energy of M identical, classical particles 
interacting pairwise: 

U M {r ll --- 1 r M ) =5>(l r i- r il)> (12) 

i<j 

where $ is the interaction potential. Now imagine a system of M + 1 particles, and consider 
two Hamiltonians: 
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M+l 2 

H A (z)= £ f- + U M {T 1 ,.--,r M ) (13) 

1=1 

M+l 2 

HB (*)= £ ^ + U M+l {v x ,.-.,v M+x ). (14) 

In both cases z is the 6(M + l)-dimensional vector specifying the positions and momenta 
of all the particles. H A describes M interacting particles, plus a "tagged" particle moving 
freely; H B describes M+l interacting particles. We take the positions to be confined 
within a cube of volume V = L 3 with periodic boundary conditions. Note that the difference 

M 

AH = H B -H A = J2 H\t M +i ~ n\) = *(ri, • • • , r M+1 ) (15) 

i=l 

is just the interaction energy between the tagged particle and the others. 

At a given temperature T, Eqs.[7] and || define partitions functions Z A and Z B , and the 
free energy difference AF, for the Hamiltonians H A and H B . This free energy difference, 
i.e. the work associated with reversibly turning on the interaction \l/ at constant V and T, 
is the excess chemical potential: 

Hex = AF (16) 

Let us now parametrize a family of Hamiltonians, H\(z), interpolating linearly from H A 
to H B : 

H X = H A + XV, (17) 

so that H = Hq and H B = H\. The intermediate Hamiltonians (0 < A < 1) represent a 
partially interacting tagged particle. 

A switching simulation then proceeds as follows. We start by sampling initial conditions 
from the canonical ensemble corresponding to H A . We then allow the (M + l)-particle 
system to evolve in time under a chosen thermostating scheme, as A is switched from to 1 
over a time r, i.e. from a non-interacting to a fully interacting tagged particle. By running 
either a single simulation with r very large (slow growth), or many simulations with smaller 
r (fast growth), we can construct an estimate of fi ex from the value (s) of external work, W T , 
performed during the simulation(s), using Eq.|| or ||. 

At the end of Section |, we alluded to the fact that, in the limit r — > 0, fast growth 
reduces to the free energy perturbation method of computing AF. For the particular case 
considered here, this is equivalent to Widom's particle insertion me thoM. which gives \i ex 
in terms of sampled values of the energy cost (AH) of suddenly materializing an additional 
particle at a randomly chosen location within the fluid. 

We performed switching simulations at various switching times r. The results of these 
will be presented in the following section. Here we describe the details and parameters of 
the simulations themselves. 

The usual Lennard- Jones interaction potential, for a pair of particles separated by a 
distance r, is given by 



$ LJ (r) = 4e 



a\ 12 /a> 6 ' 



(18) 
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where e is the depth of the potential at its minimum, and a is the van der Waals diameter 
of the particle. We instead used a modified potential: 



< r < 0.8a 

$(r) = \ $ LJ (r) + c(r - r c ) - d 0.8a <r <r c (19) 

r c < r. 




Here r c is a cutoff distance, which we took to be half the length of a side of the box: 

r c = L/2. (20) 

and the constants a, b, c, and d were chosen to preserve the continuity of $ and its derivative. 
We used $(r) rather than $Lj(r) both to avoid the singularity at r = 0, and to impose a 
strict radius of interaction equal to L/2. The need for the latter condition arose from our 
use of periodic boundary conditions, which is equivalent to an infinite lattice of identical 
cells: the requirement that $ vanish for r > L/2 guarantees that, given a pair of particles, 
each interacts only with a single image of the other. 

For the parameter values we have chosen, the difference between $(r) and $Lj(r) is very 
small for r > 0.8<r. Furthermore, we have found the core repulsion between particles to 
be sufficiently strong that two fully interacting particles will essentially never come closer 
than a distance 0.8cr of one another. These considerations suggest that the excess chemical 
potential fi ex which we plan to compute will differ little from that defined for the usual 
Lennard- Jones fluid. However, this issue is not particularly relevant for the central aim of 
this paper. 

Since the Lennard- Jones potential works especially well with inert gases, we modeled an 
argon fluid. The parameters associated with argon are@ 

a = 3.542 x 10~ 10 m = 3.542 A (21) 
e = 1.288 x 10" 21 J = 0.1854 kcal/mol (22) 
m = 6.634 x 10~ 26 kg = 39.94 amu. (23) 

We took these values to define the units of length, energy, and mass in our simulation. The 
units for time and temperature are then: a\jm/t = 2.542ps and e/ks = 93. 3K, respectively. 
Henceforth, all values will be quoted in these units. 

We chose the Andersen thermostaM to simulate the effects of an external heat reservoir. 
In this scheme, the evolution of the fluid is governed by Newton's equations, supplemented by 
the forced thermalization of the momenta of randomly chosen particles, at regular intervals. 
Specifically, at fixed time intervals St And, a particle is chosen at random, and all three 
components of its momentum vector are replaced by thermal values, sampled randomly 
from the Maxwell-Boltzmann momentum distribution at the thermostat temperature T. To 
carry out the integration of Newton's equations, we used the velocity- Verlet algorithmS. 

All of our simulations were run using M = 125 untagged particles, plus one tagged 
particle, inside a box of sides L = 5.3, with the thermostat set at T = 1. 

During each switching simulation (following a relaxation interval at A = 0, to generate 
an initial equilibrium microstate; see Section |IV|), the interaction strength of the untagged 
particle was turned on according to the schedule: 
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X(t) = (t/rf, 



(24) 



from t = to t = r. The time step for the velo city- Ver let integration was taken to be 
5ty er = .01, and the Andersen momentum thermalization was implemented at every time 
Step: St And = -01. 



III. RESULTS: ILLUSTRATION OF PRINCIPLE 

We ran a varying number of simulations of the switching process, at seven different values 
of the switching time, from N = 3334 simulations at r = 3.0, to N = 10 simulations at 
t = 1000.0, keeping the total switching time devoted to each set of simulations, Nr, fixed 
at (approximately) 10 4 . For every simulation, we computed the work performed as a result 
of turning on the interactions between the tagged particle and the rest. We then calculated 
both the ordinary and exponential average of these values of work, over the simulations 
corresponding to each value of r: 

1 N 

W?' N = E Wr, n (25) 

n=l 



W x,N = _p-i ln 



1 N 

£ eX p(-/W T , n ) 

iV n=l 



(26) 



where W T>n denotes the work performed during the n'th simulation of duration r. 

The results are summarized in Fig.|l|. The error bars on the exponential averages (circles) 
were computed with the bootstrap methodic. Error bars on the ordinary averages (squares) 
are not shown, but were typically the size of the squares themselves. 

Additionally, we performed two very long simulations at r = 25000.0, and used the 
average of these as our best estimate of the true value of AF (defined by Eq.|]). This 
is shown as the horizontal line at W = 1.174 in Fig.[I], and - judging from the difference 
between these two (nearly quasi-static) simulations - this value should be understood to 
have an uncertainty of ~ 0.1. 

Let us now consider the salient features of Fig.[TJ, beginning with the values of W® ,N . If we 
view the work performed during a single simulation, W T , to be an estimate of the free energy 
difference AF, then this estimate is subject to both statistical and systematic errorsf^rt^l. 
The former reflect the spread in the distribution of W T values (for a given r), the latter a 
bias: W T tends to over-estimate AF. This bias follows rigorously from the non-equilibrium 
work relation!, but can also be understood in terms of thermodynamic considerationsEIH 
(e.g. the Second Law), if one accepts the interpretation of W T as work. By averaging W T 
over infinitely many realizations, the statistical error is removed, but the systematic error 
remains: 



N 



lim W? N =W T > AF (finite r). (27) 



Thus, for any value of r, the quantity W T — AF represents the systematic error associated 
with using simulations of duration r to estimate AF. As expected, this error (approximated 
by W^' N — AF in Fig.[TJ) is positive, but decreases as we approach the reversible limit, r — > oo. 
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The exponential average of W T over N simulations of duration r, viewed as an estimate 
of AF, also contains errors when N is finite!, but these are removed in the limit of infinitely 
many simulations: 

lim W?' N = AF (arbitrary r). (28) 

The values of W X ' N shown in Fig.|Tj are consistent with this statement. All of the exponential 
averages closely approximate (our best estimate of) the free energy difference, AF = 1.174, 
illustrating the validity of the fast growth method. 

The histogram in Fig.[| displays the distribution of work values for the set of simulations 
performed at r = 3.0. The three vertical lines show the locations of the ordinary average 
of the work values, W a = 8.413, the exponential average, W x = 1.478, and the free energy 
difference, AF = 1.174, shown in Fig.^]. (Throughout the paper, we frequently suppress 
the subscript r and/or the superscript N.) Note that W x and AF are quite close: despite 
values of W ranging from —4.3 to +34.1, the exponential average W x differs from AF by 
only ~ 0.3, which is not much larger than our uncertainty in the value of AF itself! 

Figs.|I] and ^] are intended as an "illustration of principle" of the fast growth method, 
demonstrating that the value of AF can be extracted from an ensemble of finite-time switch- 
ing simulations, during which the system is explicitly driven away from equilibrium. 



IV. COMPARISON OF METHODS 

Both slow and fast growth converge to the correct value of AF in the limit of infinite 
computational resource: r — ► oo for slow growth, N —>■ oo for fast growth. In practice, 
of course, we can perform neither a single infinitely long simulation, nor infinitely many 
finite-time simulations, so it becomes interesting to pose the following question: Given a 
finite amount of computer time, which method is preferable? We now consider this question 
empirically, in a manner relevant to the practical computation of free energy differences. 
Specifically, we study how the accuracy of the predicted free energy difference depends on 
the choice of the switching time, r, given a fixed amount of total computational time. 

Imagine that we have been allotted ttot units of simulation time. In this total we must 
include not only the time devoted to repetitions of the switching process itself, but also 
the time used to generate a microstate sampled from equilibrium, before the start of each 
process. If each switching realization is preceded by a relaxation interval of duration trel 
(to generate the initial microstate), then the number of switching realizations possible is 
given by the largest integer N satisfying the inequality: 

N < TTOT . (29) 
trel + r 

In the simulations described in Section |TJ, for each value of r the results were obtained 
from a single, long run, in which realizations of the switching process itself alternated with 
relaxation intervals of duration trel = 1-0, during which the system evolved at fixed A = 0. 
(The final microstate at the end of a given switching simulation was taken as the seed for 
the following relaxation interval. ^1 The value of trel was chosen on the basis of simulations 
showing that correlations in our system decay to zero over ~ 0.6 time units.) We used these 
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results to investigate predictions of AF for various switching times r, given a total resource 
of ttot — 1050.0 time units of simulation. For each of the seven switching times considered, 
we generated ten predictions of AF using the data from the simulations. For instance, at 
t = 5.0, Eq.[29| gives iV = 175; therefore, of the 2000 switching realizations performed at 
t = 5.0, we used the first 1750 to generate ten independent estimates of AF. The remaining 
250 realizations were not used here. 

Fig.^ displays the results, for the various switching times. Each circle represents an 
estimate of AF obtained from a total of (no more than) 1050.0 units of simulation time. 
The data at r = 1000.0 are slow growth estimates (N = 1), whereas the remaining points 
make explicit use of the fast growth formula (Eq.|5|). A striking feature of these results is that 
the accuracy of the predictions does not seem to depend drastically on r. In other words, 
for the example considered in this paper, fast growth and slow growth yield comparably good 
estimates of AF for a fixed computational resource, over a wide range of switching times. 

Additionally, we generated ten estimates of AF using the Widom method - shown as 
crosses at r = 0.0 - again with a total of 1050.0 time units devoted to each estimate. Here, 
however, we effectively took Trel = 0.1: after every time interval At = 0.1, we generated 
a random location for the (M + l)'th particle, determined AH, and then continued with 
the simulation of M — 125 mutually interacting particles. Thus, 10500 samples of AH 
contributed to each of the ten Widom estimates. As we made no attempt to optimize At, 
these results do not really represent a direct and fair comparison between fast growth and 
the Widom method, and rather serve as an illustration of the latter. 

In the asymptotic limit of large N and large r, Crooks! has shown that the size of 
errors in the estimate of AF depends only on Nr. On the other hand, for fast switching, 
the distribution of work values can become quite wide, resulting in slow convergence of W? ,N 
with N. This is a well-known problem associated with computing exponential averages (as 
with the free energy perturbation method!!) , which manifests itself both in large statistical 
errors, and - if the distribution is wide enough - as a significant bias toward over-estimation 
of AF. The Widom data (r = 0.0) certainly shows these symptoms; the somewhat wide 
scatter of predictions at r = 3.0 is likely also evidence of this problem. 

We suspect that, generically, the expected error associated with fast growth decreases 
with increasing r (for fixed t TO t)- If this is indeed the case, then FigFJ suggests that for 
t > 5.0 this trend is lost in the noise; any gain in accuracy obtained by using one switching 
realization of duration 1000.0 rather than 175 realizations of duration 5.0, is likely to be 
smaller than the unavoidable statistical error associated with trying to compute AF with a 
finite resource of 1050.0 units of simulation time. 

The above, tentative conclusion - that for a given t TO t the accuracy does not depend 
strongly on r, over a wide range of switching times - suggests a certain advantage of fast 
growth over slow growth as a practical computational tool. Namely, fast growth allows for 
the immediate estimation of errors^, both statistical and systematic. Using the N values of 
work obtained from different switching realizations, the statistical uncertainty in the average 
of c~^ Wt - the unbiased estimate of e~ l3AF - can be ascertained in the usual way, in terms 
of the variance; the resulting upper and lower bounds defining the error bar then translate 
immediately into lower and upper bounds on the estimate of AF itself. Alternatively, one 
can use the bootstrap methode^, which is sensitive to the possibility that the exponential 
average might be dominated by one or a few particularly small values of work, out of the N 
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values obtained from the switching realizations. Slow growth, by contrast, produces only a 
single value of work, which by itself provides no estimate of its own error! 

In addition to statistical errors, the fast growth estimate contains a systematic bias: for 
finite N, W*> N tends to over-estimate AF, though to a lesser degree than VF* ,iv §. This bias 
has the same origin as the "sample-size hysteresis" studied in the context of the free energy 
perturbation method, and Eq.10 of Ref£3 gives us a leading-order estimate of its size: 



w x,n _ Ap _ p^ w /2N. (30) 

Here the overbar denotes an expectation value (i.e. the average over infinitely many attempts 
to estimate AF with N values of W, using fast growth), and a^y is the variance in the work 
values. Note that this bias vanishes as N — > oo, as demanded by consistency with Eq.[|. 

For each of the fast growth estimates shown in Fig|| we have estimated both the statis- 
tical uncertainty (using the bootstrap method) and the systematic bias (using Eq.|30|) . We 
found that the latter never exceeded 0.07, and that the estimated statistical uncertainty was 
typically about an order of magnitude greater than the estimated bias. These conclusions 
are in agreement with Fig.|3], which visually suggests that statistical rather than systematic 
errors dominate our fast growth predictions. 

With fast growth, the easy estimation of error bars can furthermore be harnessed to 
determine, "on the fly", the optimal amount of computer time to devote to the estimation 
of AF: the user simply performs one simulation after another, updating both W®' N and its 
error bar after each simulation, and stops when the size of the error becomes satisfactorily 
small. With slow growth, one must estimate ahead of time how long the simulation needs 
to be, given the desired accuracy. 

Our discussion so far has been implicitly oriented toward the computation of free energy 
differences on single-processor computers, where easy estimation of statistical errors seems 
to be the main advantage of fast growth over slow growth. However, on parallel machines 
another considerable advantage emerges: fast growth is embarrassingly parallelizable. Given 
P processors and a finite amount of dedicated run time, the user simply lets each proces- 
sor perform one or more switching simulations, and at the end all the values of work are 
gathered and the exponential average is computed. This provides the maximal efficiency 
of parallelization, ~ 100%, as there is essentially no "cross-talk" between processors. Fur- 
thermore, the additional programming effort is minimal: once a code has been written to 
perform switching simulations on a single processor, little more needs to be done to let P 
copies of the code run independently on P processors. By contrast, to efficiently distribute a 
single simulation over P processors requires considerable programming skills and time, and 
the optimal parallelization scheme may well depend on details of the hardware - e.g. are the 
processors divided into shared-memory clusters?, how fast is the communication between 
different processors?, etc. - making it less portable. 

The fast growth method is particularly well-suited for Beowulf-type clusters, built at 
relatively low cost with off-the-shelf hardware components. In such clusters, the proces- 
sors typically do not share memory, and communication between them is often relatively 
slow. Both drawbacks are avoided by any embarrassingly parallelizable method such as fast 
growth. 
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V. FAST GROWTH, SLOW GROWTH, AND LINEAR RESPONSE 



The idea of computing AF from numerous values of W T is not new. Indeed, it has been 
recognized for some time that, in the interest of having an estimate of statistical errors, 
several or more independent switching simulations should be carried out0. 

Given a number of independent values of W T , the simplest imaginable approximation of 
AF - though clearly not the best (Fig.[l|) - is just the ordinary average of these values: 

AF « W? (31) 

This amounts to slow growth, with statistical errors reduced by averaging over samples. In 
1991, Hermans^ suggested instead the following formula: 

AF^W?-{3a 2 w /2, (32) 

where a^y is the variance of the set of work values (closely related results are found in 
Wood et a@, and - for isolated systems - in Tsao et a§). This result can be viewed 
as a fluctuation-dissipation relation, and therefore ought to be valid in the linear response 
regime, i.e. when A is switched slowly enough that the system remains near equilibrium for 
the duration of the process. 

Eq.|52"| represents a correction to Eq.^T], aiming to compensate for the systematic error (see 
Section |ITT|) inherent in using W T as an estimate of AF. Thus we can expect W® — /3cr^/2 to 
give a better estimate of AF than W®, in the regime of validity of linear response. Using the 
work values obtained from the simulations described in Section |HI| , we have compared - for 
various switching times - three estimates of AF: W° (corresponding to slow growth, with the 
statistical error reduced by averaging), W® — (3a^/2 (linear response), and W* (fast growth). 
Fig.f| displays the results. As expected, there exists a near-equilibrium regime (roughly, 
t > 20.0) in which linear response does a good job of removing the systematic error inherent 
in using W T as an estimate of AF. For shorter times, however, linear response breaks down, 
as the system is driven significantly away from equilibrium. Fig.|] furthermore illustrates 
that the non-equilibrium work relation is truly a far- from- equilibrium result, remaining valid 
even outside the regime of validity of the near-equilibrium prediction (Eq.|3"2"D. 

The relationship between the three estimates considered in Fig.|] can be understood by 
rewriting Eq|J as follows!: 

A^EMr 1 ^ (33) 

n=l ll - 

where u n is the n'th cumulant of the distribution of values of W. Keeping only the first term 
on the right side gives us AF m W®, whereas keeping two terms yields the linear response 
approximation, AF « — (3 (7^/2. Therefore these approximations can be viewed as the 
first two in an expansion derived from the non-equilibrium work relation, Eq|| 

Note that one can justify keeping only the first two terms in the cumulant expansion, by 
assuming a Gaussian distribution of work values: since all higher (n > 2) cumulants vanish 
for a Gaussian, Eq.[J^ reduces to Eq.^ in that case. More directly (i.e. without resorting to 
a cumulant expansion), it is easily verified that for an assumed Gaussian distribution of work 
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values, Eq.|3] implies a relationship between the mean and the variance of the distribution, 
namely AF = W - /3a w /2. 

It is easy to argue heuristically that a Gaussian distribution of work values emerges in 
the linear response regime, i.e. when the switching is slow enough (r large enough) to keep 
the system near equilibrium. Namely, we can imagine that the work W, for one simulation, 
becomes a sum of many independent contributions: if we break the long time interval 
< t < r into many shorter segments - each, however, longer than the correlation time 
associated with the dynamics of the system - then the contributions to W T from the various 
segments will be (roughly) independent. The central limit theorem (CLT) then leads us to 
expect that, given many realizations of this process, the distribution of work values, p(W), 
will be approximately Gaussian. However, this argument requires some care: if the spread 
in the values of W is large (aw 3> ksT), then the dominant contribution to / pt~^ w comes 
from values in the (lower) tail of p(W), precisely where the CLT "breaks down"0. So a 
more careful justification invoking the CLT (to derive Eq.[32] from Eq.^) requires also that 
aw be not much larger than k^T . This will be satisfied for sufficiently slow switching, since 
p — ► 5(W — AF) as t — > oo (Eq.|I|). Hence, the argument is ultimately justified. 

CONCLUSIONS 

In this paper, we have studied the fast growth method, applying it to compute the excess 
chemical potential of a (modified) Lennard- Jones fluid. Figs.[I] and |2| illustrate the validity of 
the method, showing that it accurately estimates AF over a wide range of switching times, 
including short times for which the system is driven far from equilibrium. Moreover, Fig.|3| 
suggests that the accuracy of the prediction does not depend strongly on the choice of r (for 
a fixed total computational resource), again over a wide range of values. However, when r 
is too short, then fast growth becomes susceptible to the problem of poor convergence of 
exponential averages, familiar to any practitioner of the free energy perturbation method 
(which in the present context is the Widom particle insertion method). Finally, we have 
considered the relationship between fast growth and the linear response approximation. 

The natural question left open by our study is: How generally valid are the results 
which we have found empirically? In particular, how robust is the conclusion about relative 
performance of slow and fast growth, i.e. "comparable accuracy for equal computational 
resource"? If this conclusion is indeed generally true, then it suggests certain advantages 
for fast growth (easy estimation of errors, parallelizability) . 

Finally, it is likely that fast growth can be improved by combining it with strategies 
which have been studied in the context of the free energy perturbation (Eq.[TT|). Suggestions 
to this end have included the overlapping distributions methodl!, higher-order cumulant 
expansions!!, and averaging of forward and reverse processes!!. 
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FIGURES 



FIG. 1. For each of the seven sets of simulations, corresponding to switching times from 
r = 3.0 to 1000.0, we computed both the ordinary average (W a ) and the exponential average 
(W x ) of the work values. These are shown in the figure, along with our estimate for the true 
value of AF, the solid line at W = 1.174. The ordinary averages, W a , displayed as squares, 
over-estimate AF, although the difference decreases as one approaches the reversible limit. By 
contrast, the exponential averages W x (circles) provide significantly better estimates of AF. The 
number of simulations, N, performed for each value of the switching time, r, was chosen so that 
Nt = 10000.0 in each case (except Nt = 10002.0 for r = 3.0). 

FIG. 2. A histogram, showing the distribution of values of W, in bins of unit energy width, 
for the 3334 simulations performed at the shortest switching time, r = 3.0. The three vertical lines 
show the ordinary average of these work values, W a , the exponential average, W x , and (our best 
estimate of) the free energy difference, AF. 

FIG. 3. Each circle or cross represents an estimate of AF, using no more than 1050.0 units of 
simulation time. The ten circles at r = 1000.0 are thus slow growth estimates; all other circles are 
fast growth estimates; and the crosses at r = 0.0 were obtained with the Widom method. 

FIG. 4. The squares, circles, and diamonds show three sets of estimates of AF: W a , W x , and 
W a — /3cj^/2, respectively. These estimates were obtained from the same data used in the previous 
figures, i.e. the squares and circles are the those shown in Fig.Q. 
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Figure 1 
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Figure 3 
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